# TODO: Add comment
# 
# Author: rogb
###############################################################################

VaR <- function(Value,p=0.99){
	tmp <- quantile(Value,1-p)
	names(tmp) <- paste(p*100,"%",sep="")
	tmp
}

ES <- function(Value,p=0.99){
	tmp <- VaR(Value,p)
	sapply(tmp,function(VaR)mean(Value[Value<=VaR]))
}

EqSim <- function(mu,Cov,nProj=12,nSim=10000){
	# nSim <- 10000
	# nProj <- 12
	# mu <- 1:10/100
	# Sigma <- diag(10)/10
	nEq <- length(mu)
	simData <- array(NA,c(nEq,nProj,nSim))
	idxData <- array(NA,c(nProj+1,nEq,nSim))
	if(nEq==1){
		simData[] <- rnorm(n=nProj * nSim,mean=mu-Cov/2,sd=sqrt(Cov))
	}else{
		simData[] <- t(mvrnorm(n = nProj * nSim, mu-diag(Cov)/2, Cov)) 
	}
	idxData[] <- apply(simData,c(1,3),function(x)exp(cumsum(c(0,x))))
	idxData <- aperm(idxData,c(2,1,3))
	idxData
}

DistPlot <- function(x){
	plot(density(x))
	abline(v=ES(x))
}

SimPlot <- function(x,col=c(0,0,1),resolution=200,colScale=1,ESq=c(0.01,0.99),VaRq=c(0.01,0.99),...){
	x <- x[1,,] 
	timeMax <- dim(x)[1]-1
	time <- 0:timeMax
	
	meanSim <- apply(x,1,mean)
	if(!all(is.na(ESq))){
		esSim <- matrix(NA,ncol=(timeMax+1),nrow=length(ESq))
		esSim[] <- apply(x,1,function(x,ESq){
					if(all(ESq>=0.5)){
						out <- ES(x,ESq)
					}
					if(all(ESq<0.5)){
						out <- -ES(-x,1-ESq)
					}
					if(any(ESq<0.5) & any(ESq>=0.5)){
						out <- c(-ES(-x,1-ESq[ESq<0.5]),ES(x,ESq[ESq>=0.5]))
					}
					out
				},ESq)
	}
	
	if(!all(is.na(VaRq))){
		varSim <- matrix(NA,ncol=(timeMax+1),nrow=length(VaRq))
		varSim[] <- apply(x,1,quantile,VaRq)
	}

			
	nSteps <- floor(resolution/2)*2
	qSim <- t(apply(x,1,quantile,seq(0,1,length.out=nSteps)))
	
	plot(NULL,xlim=c(0,timeMax),ylim=range(x),...)
	iMax <- (nSteps) / 2 
	for(i in 1:iMax){
		# i = 1
		polygon(c(time,rev(time)),c(qSim[,i],rev(qSim[,i+1])),col=rgb(col[1],col[2],col[3],(i/iMax)^colScale),border = NA)
		polygon(c(time,rev(time)),c(qSim[,nSteps+1-i],rev(qSim[,nSteps-i])),col=rgb(col[1],col[2],col[3],(i/iMax)^colScale),border = NA)
	}	
	grid(NULL,NULL)
	lines(time,meanSim,col=1,lwd=2,lty=2)
	if(!all(is.na(ESq))){
		for(i in 1:dim(esSim)[1]){
			lines(time,esSim[i,],col=2,lwd=2,lty=2)
		}
	}
	if(!all(is.na(VaRq))){
		for(i in 1:dim(varSim)[1]){
			lines(time,varSim[i,],col=3,lwd=2,lty=2)
		}
	}
	if(all(is.na(VaRq)) & all(is.na(ESq))){
		legendText <- c("Mean")
		legendCol <- c(1)
	}
	if(all(is.na(VaRq)) & all(!is.na(ESq))){
		legendText <- c("Mean",paste(paste(ESq*100,"%",sep="",collapse=",")," Expected Shortfall",sep=""))
		legendCol <- c(1,2)
	}
	if(all(!is.na(VaRq)) & all(is.na(ESq))){
		legendText <- c("Mean",paste(paste(VaRq*100,"%",sep="",collapse=",")," Value at Risk",sep=""))
		legendCol <- c(1,3)
	}
	if(all(!is.na(VaRq)) & all(!is.na(ESq))){
		legendText <- c("Mean",paste(paste(ESq*100,"%",sep="",collapse=",")," Expected Shortfall",sep=""),paste(paste(VaRq*100,"%",sep="",collapse=",")," Value at Risk",sep=""))
		legendCol <- c(1,2,3)
	}
	legend("topleft",legendText,col=legendCol,lty=2,lwd=2)
}

